Fast Diffusion Mechanism of Silicon Tri-interstitial Defects 
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Molecular dynamics combined with the nudged elastic band method reveals the microscopic self- 
diffusion process of compact silicon tri-interstitials. Tight-binding molecular dynamics paired with 
ab initio density functional calculations speed the identification of diffusion mechanisms. The dif- 
fusion pathway can be visualized as a five defect-atom object both translating and rotating in a 
screw- like motion along (111) directions. Density functional theory yields a diffusion constant of 
4 x 10 _o exp(— 0.49 eV/fesT) cm 2 /s. The low-diffusion barrier of the compact tri-interstitial may 
be important in the growth of ion-implantation-induced extended interstitial defects. 
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Following high-energy ion implantation, strongly su- 
persaturated silicon self-interstitials can agglomerate to 
form macroscopic planar interstitial structures, {311} de- 
fects pj . High-temperature annealing is necessary to re- 
move lattice damage following ion implantation. How- 
ever during the annealing process the {311} defects com- 
prise a large reservoir of interstitials 0, 0, 3 , which are 
released upon annealing and hence drive boron transient 
enhanced diffusion 0,00' an undesirable process which 
causes spatial broadening of boron concentration profiles. 
On the other hand, following low-energy implantation, 
significant boron TED is still observed, even though no 
visible {311} defects are developed Q- 

It is crucial to understand the diffusion of various 
point defects in order to have a quantitative under- 
standing of boron TED. Measurements of the diffusion 
rate of defects in silicon have been reported in experi- 
ments H H EIH |H El. Current experimental 
techniques cannot cleanly resolve the diffusion rates of 
more complex defect species, and due to the atomic-scale 
size of point defects, defect diffusion pathways cannot be 
resolved at all [l5j . Thus, numerical simulations provide 
a unique way to study the technologically-important dy- 
namics of point defects in crystalline silicon. 

Previously, such simulations have studied numerous 
defect species u sing both classical and quantum Hamil- 
tonians 0, 0, llSL [Hi l2(i| . Using nudged elastic band 
methods (NEB) [21j within density functional theory 
(DFT) code, Lopez et al. compute an activation en- 
ergy of 0.28 eV for a neutral single interstitial diffusing 
along the X- T-X path (2^ • Kim et al. estimates the self- 
diffusion barrier of a particular di-interstitial performing 
a reorientation to be 0.5 eV 17]. Recent DFT-NEB cal- 
culations done in our group confirm this pathway but 
give a somewhat lower barrier of 0.3 eV |38j . 

Cogoni et al. 23] use temperature-accelerated 
molecular-dynamics [24[ with the Kwon et al. 2f| tight- 
binding (TB) potential to systematically study the dif- 
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FIG. 1: Diffusion of if to the neighboring site along the 
(111) direction. The initial configurations viewed from [110] 
and [111] are shown in (a) and (b), respectively. The final 
configurations viewed from [110] and [111] are shown in (c) 
and (d), respectively. Defect atoms are shaded. All of the 
above structures are relaxed using density functional theory. 



fusion of low- lying single-, di- and tri-interstitials. In 
particular they find all are fast diffusers with diffusion 
barriers of 0.94, 0.89 and 1.71 eV, respectively. We ex- 
tend that work both with a more accurate tight-binding 
potential and with density-functional theory to further 
refine the diffusion pathway and the diffusion constant. 
Here we concentrate on the tri-interstitial finding an ac- 
tivation energy of 0.4-0.5 eV in DFT (and 0.6 eV in TB). 

Which tri-interstitial to study? In a previous work 
from our research group |26| . the three lowest-energy 
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Reaction Coordinate 

FIG. 2: Transition paths of it in DFT-GGA and the approx- 
imation of collective motion of five atoms within DFT-GGA. 
Five atoms are highlighted in gray for the initial state, saddle 
point and final state in the insets, respectively. 



tri-interstitial geometries, denoted 1$, i|, and I§, were 
identified by analysis of tight-binding molecular dynam- 
ics (MD) simulations followed by density functional re- 
laxations. The TB calculations use the Lenosky et al. [27J 
potential. If is the ground state and /f and If are 
excited states j3^; the density functional energies for 
the three defects in 216+3 atom supercells were 2.24 
eV/atom, 2.37 eV/atom, and 2.47 eV/atom respectively. 
The tight-binding MD simulations show only 7| is a rapid 
diffuscr, while 1^ and J| can be formed by interconver- 
sion of 7g, but are themselves immobile within the 5 ns 
time scale of the simulations |26j. Hence, we focus on 
elucidating the microscopic diffusion process of the com- 
pact tri-interstitial J|. i| has a compact defect geometry 
in which a perfect tetrahedron of four atoms replaces a 
single atom in the silicon lattice, with the faces oriented 
along the four symmetry related (111) directions. 

Tight-binding molecular dynamics determines an ini- 
tial /| diffusion path. Density-functional theory uses 
the nudged-elastic-band method to refine the diffusion 
path and determines the accurate diffusion barrier. The 
self-diffusion of j| can be visualized in terms of the 
atoms most distorted during the process. In partic- 
ular a five-atom object translates while rotating to 
avoid adjacent atoms. NEB finds four equivalent paths 
along (111) with a diffusion constant of D = 4 x 
1(T 5 exp(-0.49 eV/k B T) cm 2 /s. 

Minimum energy path / nudged elastic band. Figure ^ 
shows two /| defects in neighboring sites that correspond 
to initial and final configurations along the diffusion path, 
the [111] direction. We obtain the initial pathway from 
analysis of molecular dynamics trajectories 26] . The 
climbing-image NEB (CI-NEB) method ^ refines an 
accurate pathway and finds the diffusion barrier between 
these two i| minima. The CI-NEB scheme guarantees 
that the image with the highest energy converges to the 
saddle point. All NEB calculations are performed with 




FIG. 3: Five most displaced atoms during the diffusion. 
These five atoms move collectively. Atom D and E define 
an axis of rotation ([111]), around which atom A, B and C 
are rotating. 

the VASP density functional code E^, |^3] employing 
the Perdew-Wang GGA functional J3ll|. and ultra-soft 
Vanderbilt-type pseudo-potentials [33 ] as provided by G. 
Kresse and J. Hafner [33|]. An initial relaxation of the 
two end points initiates a full relaxation of seven im- 
ages along the path, keeping the volume of the cell fixed. 
Energy and atomic position convergence of 3 meV and 
0.005 A, respectively, is confirmed for 64+3 atom super- 
cell by comparing the results for a 250 eV energy cutoff 
and a 3 x 3 x 3 k-point mesh with a 300 eV cutoff and 
a 4 x 4 x 4 mesh. The seven images are initialized by 
linear interpolation between the two relaxed end points. 
Each of the images is relaxed until the atomic and spring 
forces are less than 10 me v/A. 

Figure [5] shows the diffusion path with an activation 
energy of 0.49 eV. Harmonic transition state theory |34[ 
calculates the defect jump rate T from the phonon fre- 
quencies at the minimum z/™ m and at the saddle point 
vf ad and the activation energy AE 

r = r e - AE ' ksT , (l) 

where the prefactor Tq is given by: 

3N-3 3N-4 

r = n < in i n v t ad - (2) 

i i 

The dynamical matrix method within ab initio GGA den- 
sity functional theory determines the phonon frequen- 
cies. Each atom is displaced in the x, y, z direction by 
0.03 A and the calculated forces are used to construct 
the Hessian matrix for the system, which is diagonalized 
to find the phonon frequencies. This yields a jump rate 
T = 0.2 THz exp(-0.49/fc B T). Calculations for a larger 
cell of 216 + 3 atoms provide a diffusion barrier of 0.43 eV 
estimating a finite size error of about 0.1 eV. 

Collective motion of atoms. In the diffusion event, the 
five most-displaced atoms move collectively, with a screw- 
like motion. The insets in Figure |21 showing the atomic 
configurations for initial, saddle and final states, indicate 
significant displacement of the five solid-colored atoms. 
All other atoms relax only slightly during the diffusion 
with a maximum displacement of 0.18 A. 

In order to illustrate the screw-like collective motion, 
we highlight the five highly-displaced atoms in Figure 
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Translation (in unit of 1.31 A) 

FIG. 4: Transition path with seven images projected onto 
two-dimensional reduced coordinates. The inset shows images 
3, 5 (saddle point), and 7 along the transition path. 



TABLE I: The translation and rotation of five atoms along 
the [111] direction with respect to the initial state for all seven 
images in Fig 2] 



image 


Translation 


(A) Rotation (degree) 


1 


0.18 


0.4 


2 


0.42 


4.6 


3 


0.57 


16.4 


4 


0.66 


30.0 


5 


0.75 


43.6 


6 


0.89 


55.4 


7 


1.13 


59.6 



Atom A, B, and C are located in a (111) plane, forming 
an equilateral triangle, and translate along and rotate 
about the [111] axis defined by atom D and E. Mean- 
while, atoms D and E translate along the same [111] 
direction. Atoms D, A and B form an equilateral tri- 
angle, whose normal vector indicates another diffusion 
direction, and all five atoms form a double-tetrahedron. 
The bond between A and B is 2.49 A at the minimum, 
6% longer than the perfect silicon-silicon bulk bond dis- 
tance. During the diffusion event this bond length varies 
by less than 4.5%. During defect migration the double- 
tetrahedron translates 1.31 A and rotates 60 degrees as 
shown in Figm-e^and the insets of Figure^ TableHJlists 
the values for the translation and rotation of the seven 
relaxed NEB images. 

The diffusion event is approximated by a uniform 
translation and rotation of a fixed double-tetrahedron. 
For instance, rotating the double-tetrahedron of the ini- 
tial configuration (image 1) by 16.4 degrees, and trans- 
lating it by 0.57 A while keeping its geometry fixed and 
relaxing the other atoms provides an approximation to 
image 3. Applying this procedure to each of the im- 



TABLE II: Formation and activation energy of the neutral 
single-interstitial diffusing along the X-T-X path, and of the 
neutral ground state di-interstitial and compact tri-interstitial 
during self-diffusion. The formation energies are calculated 
in a 216 + n atom supercell and the activation energies are 
calculated using NEB with a 64 + n atom supercell. Details 
of the diffusion paths of single- and di-interstitials will be 
published elsewhere [38|. 

Formation energy /atom (eV) Activation energy (eV) 
7i 3^10 0.29, 0.28 a 
h 2.827 0.30 
_Z| 2.368 049 

"NEB result from Ref. |22| is calculated using 216 + 1 atom su- 
percell. 



ages along the diffusion path, we obtain the approximate 
diffusion path in Figure [5] The approximate activation 
energy is 0.6 eV, only 20% higher than the fully relaxed 
value. 

Figure 01 shows the translation and rotation of the 
double-tetrahedron during the diffusion. Translation 
occurs before the structure rotates. The double- 
tetrahedron is displaced half way and rotated 30 degrees 
at the saddle point. The insets show images 3, 5 (saddle 
point) and 7 viewed from the [111] direction and illustrate 
that three-fold symmetry persists during the diffusion. 

Diffusion rate. The defect J| can move to four neigh- 
boring sites along the four symmetry-equivalent (111) di- 
rections. The diffusion constant for this random walk is 
D = 2 |a 2 r (derivation is given in Ref. |35j]h where 
the jump rate T is defined by Eq. ^ and the displace- 
ment during the diffusion event is a — 1.31 A. The fac- 
tor two takes into account the fact that the defect can 
rotate either in clockwise or counter-clockwise direction 
during the diffusion. The resulting diffusion constant is 
D = Ax 10- 5 exp(-0.49 eV/k B T) cm 2 /s. 

Discussion. Table [H] summarizes the ab initio den- 
sity functional theory results for formation and activa- 
tion energy of I\, I2, and compact ^3 defects. The the- 
oretical values of 0.28 to 0.49 eV are considerably lower 
than the experimental results from metal-diffusion ex- 
periments which find activation energies, believed to be 
for single intcrstitials, typically ranging from ~1.4 to 
1.8 eV j3||. However, a recent experiment 0] identi- 
fies enhanced diffusion of single interstitials at 150 K, 
while noting it is inconsistent with results from earlier 
approaches. This lower temperature would be consistent 
with the NEB result for the energy barrier of single in- 
terstitial diffusion of 0.3 eV. 

Our pathway for ^3 diffusion is the same as that found 
by Cogoni et al. [23] in the 64 + 3 atom supercell. On 
the other hand, our DFT barriers for I\, I2, and J3 of 
0.3, 0.3 and 0.4-0.5 eV are much smaller than the tight- 
binding results by Cogoni et al. 23] of 0.9, 0.9 and 1.7 
eV. We also perform the NEB calculation to examine 
the diffusion path of the compact tri-interstitial within 
Kwon's potential [2{|, and obtain a diffusion barrier of 
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1.2 eV in contrast to 0.6 eV barrier within Lenosky's TB 
potential [27|]. The Kwon's potential overestimates the 
phonon frequencies |2*t| . which makes a local minimum 
steeper. This suggests that Kwon's potential will conse- 
quently overestimate the barrier for the interstitial diffu- 
sion, which is characterized by the small displacements 
of several atoms deviating from the equilibrium sites. 

We perform density of states (DOS) calculations for 
the interstitial defect along the diffusion path within a 
216 + 3 atom supercell. The DOS is calculated for the 
minimum, the saddle point and an intermediate struc- 
ture. The DOS shows a band gap of 0.71-0.74 eV for all 
three configurations along the diffusion pathway. There 
are no states in the gap. The size of the gap is close to the 
band gap of pure Si in GGA of 0.72 eV. The lack of gap 
states indicates that charge transfer may not play a sig- 
nificant role for the diffusion of a compact tri-interstitial. 

While the formation energy of -Z3 is lower than I\ on a 
per-atom basis, it is much higher when considered on a 
total basis. One implication that can be drawn from this 
is that under equilibrium conditions, 7 3 is not present 
in any significant quantity, but under conditions which 



inject excess interstitials 5j, such as ion implantation, 
it may be present in significant amounts. Our results 
have possible relevance for modeling transient enhanced 
diffusion which follows ion implantation |37| . 

Conclusions. We have elucidated the pathway for 
the diffusion of the compact tri-interstitial, 7|, the only 
fast-diffusing tri-interstitial species in crystalline silicon. 
During the self-diffusion event, five atoms move collec- 
tively in a screw-like motion along one of four symmetry- 
related (111) directions. Our DFT result shows a low 
activation energy of 0.49 eV and a diffusion constant of 
4 x 10~ 5 exp(-0.49 eV/k B T) cm 2 /s. Under conditions 
such as ion implantation that creates excess interstitials 
and hence favor cluster formation, i| diffusion may be 
an important process due to the low activation energy. 
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